The following files will be used in this lab, all available on Canvas:

  • Climate dataset for Western North America: WNAclimate.csv
  • New York state polygon data in a shapefile: NY_Data.zip
  • A NetCDF file of global monthly air temperature: air.mon.ltm.nc

You will need to download these files from Canvas, and move them from your Downloads folder to the datafiles folder that you made previously. Make sure to unzip the zip files so that R can access the content. Note that on Windows, you will need to right-click on the file and select ‘Extract files’.

You will also need to make sure the following packages are installed on your computer:

pkgs <- c("mapview",
          "raster",
          "RColorBrewer"
          "sf",
          "tmap",
          "viridis")

install.packages(pkgs)
library(mapview)
library(raster)
library(RColorBrewer)
library(sf)
library(tmap)
library(viridis)

Static Maps with tmap

tmap is built on top of ggplot2 and works in a similar way by building a series of layers with map geometries and elements. We start by using tm_shape() to identify the spatial object to be used, and then geometries are added, including filled polygons, borders, legends, etc.


Polygon data

We’ll start by making some maps based on the NY shapefile. Start by loading this, then extract only the polygons belonging to Syracuse.

NY8 <- st_read("../datafiles/NY_data/NY8_utm18.shp")
## Reading layer `NY8_utm18' from data source 
##   `/Users/u0784726/Dropbox/Data/devtools/geog6000/datafiles/NY_data/NY8_utm18.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 281 features and 17 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 358241.9 ymin: 4649755 xmax: 480393.1 ymax: 4808545
## Projected CRS: WGS 84 / UTM zone 18N
Syracuse <- NY8[NY8$AREANAME == "Syracuse city", ]

First, let’s make a simple map showing the polygon outlines using tm_borders():

tm_shape(Syracuse) + tm_borders()

The function tm_fill() can be used to define the polygon fill color based on an attribute in the Syracuse data set (POP8). Note that this automatically adds a legend within the frame of the figure:

tm_shape(Syracuse) + 
  tm_borders() + 
  tm_fill("POP8")

The color scale can be changed by setting the palette argument in tm_fill(). This includes the ColorBrewer scales described above,and the different intervals. For example, to use the ‘Greens’ palette with percentile breaks:

tm_shape(Syracuse) + 
  tm_borders() + 
  tm_fill("POP8", palette = "viridis", style = "quantile")

Other map elements can be added. Here we add a longitude/latitude graticule with tm_graticules(), a north arrow with tm_compass(), and a line of text showing the date the map was made with tm_credits().

tm_shape(Syracuse) + 
  tm_graticules(col = "lightgray") + 
  tm_borders() + 
  tm_fill("POP8", palette = "Greens", style = "quantile") + 
  tm_compass(position = c("left", "bottom")) + 
  tm_credits("2019-10-19", position = c("right", "top"))

Point data

We’ll next make some maps with the Western North American site data. As before, we load the data, then convert to an sf object. We also load the shapefile of country outlines:

wna_climate <- read.csv("../datafiles/WNAclimate.csv")

wna_climate <- st_as_sf(wna_climate, 
                        coords = c("LONDD", "LATDD"),
                        crs = 4326)

countries <- st_read("../datafiles/ne_50m_admin_0_countries/ne_50m_admin_0_countries.shp", 
                     quiet = TRUE)

Individual symbols can be plotted on a color scale using tm_symbols.

tm_shape(wna_climate) + tm_symbols(col = "Jan_Tmp")

This takes the same arguments as tm_fill() for the color palette. We’ll use a red to blue color scale from RColorBrewer. The minus sign before the palette name reverses the order of the colors. As there is a large amount of overlap between the sites, we also add an alpha level to make the symbols transparent.

tm_shape(wna_climate) + 
  tm_symbols(col = "Jan_Tmp", alpha = 0.5, palette = "-RdBu")

We’ll next add country boundaries from the Natural Earth shapefile loaded earlier. Note that as this is a different spatial object, we have to use tm_shape() a second time to reference this, then use tm_borders() to add the lines.

tm_shape(wna_climate) + 
  tm_symbols(col = "Jan_Tmp", alpha = 0.75, palette = "-RdBu") + 
  tm_shape(countries) + 
  tm_borders(col = "gray")

Note that here we’ve plotted the symbols first to set the plotting window. Unfortunately, this means that the country lines are added on top of the symbols. If we add the country outlines first, then the window is set to the globe.

tm_shape(countries) + 
  tm_borders(col = "gray") +
  tm_shape(wna_climate) + 
  tm_symbols(col = "Jan_Tmp", alpha = 0.75, palette = "-RdBu")

To get around this, we can set the bounding box for the countries object using the bounding box of the point data. We’ll also move the legend to a better position in the map, and add a title

tm_shape(countries, bbox = st_bbox(wna_climate)) + 
  tm_borders(col = "gray") +
  tm_shape(wna_climate) + 
  tm_symbols(col = "Jan_Tmp", alpha = 0.75, palette = "-RdBu") +
  tm_layout(main.title = "W. North America Climate", 
            legend.position = c("left", "bottom"))

We can also use one of tmap’s preset themes with tm_style():

tm_shape(countries, bbox = st_bbox(wna_climate)) + 
  tm_style("cobalt") + 
  tm_borders(col = "gray") +
  tm_shape(wna_climate) + 
  tm_symbols(col = "Jan_Tmp", alpha = 0.75, palette = "-RdBu") +
  tm_layout(main.title = "W. North America Climate", 
            main.title.size = 1, 
            legend.position = c("left", "bottom"))

Finally, we can write this map to a file using tmap_save(). To do this, assign the output of all the tmap functions to an R object. Then specify this object and a file name. I’ve used pdf output here, but you can also use png, tif or most other common format extensions.

map1 <- tm_shape(countries, bbox = st_bbox(wna_climate)) + 
  tm_borders(col = "gray") +
  tm_shape(wna_climate) + 
  tm_symbols(col = "Jan_Tmp", alpha = 0.75, palette = "-RdBu") +
  tm_layout(main.title = "W. North America Climate", 
            legend.position = c("left", "bottom")) + 
  tm_style("cobalt")

tmap_save(wna_map, "wna_jan_tmp.pdf")

Facetted maps

tmap includes a similar facetting function to ggplot2. This allows you to make several different versions of the same map split by some grouping factor. A good example of this is where you have repeated measures over time for a set of locations, and you would like to display each time periods separately. You can make all the individual maps, and then add them together with tmap_arrange(), or use tm_facets(). To illustrate this, we’ll load a built in dataset (metro) that has population sizes for major metropolitan centers

data(metro)
head(metro)
## Simple feature collection with 6 features and 12 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -64.18105 ymin: -34.60508 xmax: 69.17246 ymax: 36.7525
## Geodetic CRS:  WGS 84
##            name             name_long iso_a3 pop1950 pop1960 pop1970 pop1980
## 2         Kabul                 Kabul    AFG  170784  285352  471891  977824
## 8       Algiers El Djazair  (Algiers)    DZA  516450  871636 1281127 1621442
## 13       Luanda                Luanda    AGO  138413  219427  459225  771349
## 16 Buenos Aires          Buenos Aires    ARG 5097612 6597634 8104621 9422362
## 17      Cordoba               Cordoba    ARG  429249  605309  809794 1009521
## 25      Rosario               Rosario    ARG  554483  671349  816230  953491
##     pop1990  pop2000  pop2010  pop2020  pop2030                    geometry
## 2   1549320  2401109  3722320  5721697  8279607   POINT (69.17246 34.52889)
## 8   1797068  2140577  2432023  2835218  3404575     POINT (3.04197 36.7525)
## 13  1390240  2591388  4508434  6836849 10428756   POINT (13.23432 -8.83682)
## 16 10513284 12406780 14245871 15894307 16956491 POINT (-58.40037 -34.60508)
## 17  1200168  1347561  1459268  1562509  1718192  POINT (-64.18105 -31.4135)
## 25  1083819  1152387  1298073  1453814  1606993 POINT (-60.63932 -32.94682)

This is an sf object, with a set of columns with population size for each decade between 1950 and 2030. To plot this as a facetted map, we first need to convert this to long format with one column for area, one for population size and one for year. We’ll do this imply by concatenating the different populations sizes, and repeating the location names and year, then put it all into a temporary data frame:

pop <- c(metro$pop1950, metro$pop1960, metro$pop1970,
         metro$pop1980, metro$pop1990, metro$pop2000)

name <- rep(metro$name, 6)
year <- rep(seq(1950, 2000, by = 10), each = nrow(metro))

tmp <- data.frame(pop, name, year)

Now we’ll merge this with the original sf object using the location names:

metro_long <- merge(metro, tmp, by = "name")

The merge() function acts in a similar way to SQL’s JOIN function. The name= argument allows you to specify which columns are to be used to merge.

Now we can go ahead and plot this. We’ll use symbol size to represent the population size and facet by year:

tm_shape(metro_long) +
  tm_facets(by = "year") +
  tm_symbols("pop", col = "darkred") +
  tm_shape(countries) +
  tm_borders() +
  tm_layout(main.title = "Metropolitan center population")


Raster data

In the final example, we’ll make figures using the global air temperature dataset. Start by re-reading the data (we also change the name of the layer, to reduce the amount of typing later on…).

air_temp <- rotate(raster("../datafiles/air.mon.ltm.nc", varname = "air"))

air_temp
## class      : RasterLayer 
## dimensions : 73, 144, 10512  (nrow, ncol, ncell)
## resolution : 2.5, 2.5  (x, y)
## extent     : -178.75, 181.25, -91.25, 91.25  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : memory
## names      : Monthly.Long.Term.Mean.Air.Temperature.at.sigma.level.0.995 
## values     : -37.77505, 33.58953  (min, max)
## z-value    : 0000-12-30
names(air_temp) <- "jan_tmp"

proj4string(air_temp) <- CRS("+init=epsg:4326")

Now let’s make a plot using tm_raster(), which again takes similar options for color palettes. We’ll also add the country borders.

tm_shape(air_temp) + 
  tm_raster(col = "jan_tmp", style = "fisher", palette = "-RdBu") +
  tm_shape(countries) + 
  tm_borders()

We can improve this a little by moving the color legend outside of the plotting area. We’ll increase the number of color classes to 9, and add a histogram showing the frequency of different values.

tm_shape(air_temp) + 
  tm_raster("jan_tmp", 
            style = "fisher", 
            palette = "-RdBu", 
            legend.hist = TRUE, 
            n = 9) +
  tm_shape(countries) + 
  tm_borders() + 
  tm_layout(legend.outside = TRUE, legend.outside.position = "left")


Tips and Tricks

For helpful tips and tricks, try:

tmap_tip()
## The filter argument of tm_shape can be used to filter spatial units. 
## New since tmap 2.0 
## 
## data(World)
## tm_shape(World, filter = World$continent == "Africa") +
## tm_polygons("HPI")


Interactive maps with mapview

The mapview package provides functions to produce interactive maps in R by building on the R leaflet package, which provides a low-level interface to the javascript library of the same name. mapview supports multiple spatial data types, including points, lines, polygons, and rasters.

mapview(wna_climate, 
        col.regions = "darkred", # fill color
        color = "gray",          # outline color
        alpha.regions = 0.2,     # fill transparency
        alpha = 0.3)             # outline transparency
mapview(Syracuse, 
        col.regions = "darkgreen",
        color = "black")
fn <- system.file("extdata", "kiliNDVI.tif", package = "mapview")

kilimanjaro_NDVI <- raster::stack(fn)[[1]]

mapview(kilimanjaro_NDVI,
        col.regions = viridis::viridis(n = 10)) # supply custom color palette


Here we will only review a few of the many features that mapview provides.

Basemaps

The default basemap called when using mapview() is CartoDB.positron. To use other third-party map tiles, we use the map.types argument in mapview(). For a list of basemap services, see here. For convenience, the leaflet package also provides a named vector of these tile services supported by the leaflet plugin. The vector is called providers and you can access elements of the vector like you would any other:

leaflet::providers[[35]]
## [1] "MapBox"
leaflet::providers$CartoDB.Positron # for the minimalists out there
## [1] "CartoDB.Positron"
mapview(Syracuse, 
        col.regions = "darkgreen",
        color = "black",
        map.types = c("Stamen.Watercolor", "Stamen.TonerLite"))


Layers

mapview provides for layering much in the same way that ggplot2 does.

Syracuse_centers <- st_centroid(Syracuse) 

lyr1 <- mapview(Syracuse, 
                col.regions = "darkgreen",
                color = "black")

lyr2 <- mapview(Syracuse_centers,
                col.regions = "darkblue",
                cex = 3) # size of point

lyr1 + lyr2


Aesthetic Mapping

We can also map attributes of the spatial data to various aesthetics.

mapview(Syracuse,
        zcol = "POP8",
        col.regions = viridis::cividis(n = 10),
        alpha.regions = 0.85)
# breweries in Franconia, Germany
franconia_breweries <- mapview::breweries

# set point size to number of different types of beer served at each brewery
mapview(franconia_breweries, cex = "number.of.types")


Inset Map

For whatever reason, mapview does not currently support interactive inset maps. However, these can still be added using leaflet and a little tinkering with the mapview object.

add_miniMap <- function(x) leaflet::addMiniMap(x@map) 

imap <- mapview(franconia_breweries, cex = "number.of.types")

add_miniMap(imap)